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Abstract 

We show how to solve hyperbolic equations numerically on unbounded domains by compactifi- 
cation, thereby avoiding the introduction of an artificial outer boundary. The essential ingredient 
is a suitable transformation of the time coordinate in combination with spatial compactification. 
We construct a new layer method based on this idea, called the hyperboloidal layer. The method 
is demonstrated on numerical tests including the one dimensional Maxwell equations using fi- 
nite differences and the three dimensional wave equation with and without nonlinear source terms 
using spectral techniques. 

Keywords: Transparent (nonreflecting, absorbing) boundary conditions, perfectly matched 
layers, hyperboloidal layers, hyperboloidal compactification, wave equations, Maxwell 
equations. 



1. Introduction 

Hyperbolic equations typically admit wavelike solutions that oscillate infinitely many times 
in an unbounded domain. Take a plane wave in one spatial dimension with frequency a> and 
wave number k, 

U(x,t) = (1) 

Any mapping of such an oscillatory solution from an infinite domain to a finite domain results in 
infinitely many oscillations near the domain boundary, which can not be resolved numerically. 
We refer to this phenomenon as the compactification problem JT]. It is commonly stated that 
hyperbolic partial differential equations are not compatible with compactification, and therefore 
can not be solved on unbounded domains accurately. 

A suitable transformation of the time coordinate, however, leads to a finite number of oscil- 
lations in an infinite spatial domain. Introduce 



t(x, t) - t - 




where C is a finite, positive constant. The plane wave ([TJ becomes 

u(x,t) = g-2»'X*c/(i+x)+^) _ (3) 

This representation of the plane wave has only k C cycles along a constant time hypersurface in 
the unbounded space x € [0, oo), and is therefore compatible with compactification. 
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The simple idea just described has far reaching consequences. In numerical calculations of 
hyperbolic equations one typically truncates the unbounded solution domain by introducing an 
artificial outer boundary that is not part of the original problem. Boundary conditions, called 
transparent, absorbing, radiative, or nonreflecting, are constructed to simulate transparency of 
this artificial outer boundary. There has been significant developments in the treatment of ar- 
tificial outer boundaries since the 70s, but there is no consensus on an optimal method 12] [3) . 
Especially the construction of boundary conditions for nonlinear problems is difficult (4). A suc- 
cessful technique for numerical calculations on unbounded domains resolves this problem for 
suitable hyperbolic equations and provides direct quantitative access to asymptotic properties of 
solutions. 

Furthermore, the numerical construction of oscillatory solutions as Q can be very efficient. 
Numerical accuracy requirements for hyperbolic equations are typically given in terms of num- 
bers of grid points per wavelength. In the example presented above, the free parameter C de- 
termines the number of cycles to be resolved, which can be chosen small. This suggests that 
high order numerical discretizations requiring a few points per wavelength can be very efficient 
in combination with time transformations of the type (|2ji. 

The rest of the paper is devoted to the discussion of time transformation and compactifica- 
tion for hyperbolic equations. The theoretical part of the paper (sections |2j and [3]) includes a 



detailed description of the method. We discuss the compactification problem (section 2.1 1 and 



its resolution (section |2T2| ) for the advection equation in one dimension. In section 2.3 we discuss 



the wave equation with incoming and outgoing characteristics. We show that the method works 



also for systems of equations (section 2.4 1. Hyperboloidal layers are introduced in section 2.5 
in analogy to absorbing layers. In multiple spatial dimensions, compactification is performed 
in the outgoing direction in combination with rescaling to take care of the asymptotic behavior 



(sections 3.1 and 3.2 1. The layer strategy in multiple dimensions allows us to employ arbitrary 
coordinates in an inner domain, where sources or scatterers with irregular shapes may be present 
(section [33). We finish the theoretical part discussing possible generalizations of the method to 



nonspherical coordinate systems (section 3.4 1. Section|4]includes numerical experiments in one 



and three spatial dimensions. In one dimension, we solve the Maxwell equations using finite 



difference methods (section 4.1 1. A stringent test of the method is the evolution of off-centered 
initial data for the wave equation in three spatial dimensions with and without nonlinear source 
terms (section |4~2] >. We conclude with a discussion and an outlook in section|5] 

2. Compactification in one spatial dimension 

2.1. Spatial compactification 

Consider the initial boundary value problem for the advection equation 

d t u + d x u = 0, u(x,0) - u Q (x), u(0,t) - b(t). (4) 

The problem is posed on the unbounded domain x e [0, oo). We transform the infinite domain in 
x to a finite domain by introducing the compactifying coordinate p via 

P(x) = x{p) = (5) 

1 + x 1 - p 

The advection equation becomes 

d t u + {\-p) 2 d p u = 0. (6) 
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Figure 1: Characteristic diagram for the advection equation after spatial compactification. The characteristic speed 
approaches zero near spatial infinity at jp = 1 } causing loss of numerical resolution: this is the compactification problem. 

The spatial domain is now given by p e [0,1]. Characteristics of this equation are solutions to 
the ordinary differential equation 



dp(t) 
dt 



-{I- pit)) 2 



They are plotted in figure[T] The compactification problem is clearly visible: the coordinate speed 
of characteristics approaches zero near a neighborhood of the point that corresponds to spatial 
infinity. The advection equation has a finite speed of propagation, therefore its characteristics 
can not reach infinity in finite time. 

A concrete example illustrates the problem for oscillatory solutions. Set initial data uq(x) = 
sin(27rx) and boundary data b(t) = - sin(27rf) in Q. We obtain the solution 



u(x, t) = sin(27r(^; - tj), 
which reads in the compactifying coordinate |5]) 

P 



u(p, t) = sin 2n 



1-p 



- 1 



(7) 



(8) 



The solution is depicted in figure [2] at t — on x e [0, 10] in the original coordinate and on 



1.0 

0.5 

-0.5 
-1.0 




Figure 2: The solution at time t = is plotted on the left panel. The same solution in the compactifying coordinate 
given in JSJ plotted on the right panel illustrates infinite blueshift in frequency. 

p G [0, 10/11] in the compactifying coordinate. The oscillations can not be resolved in the 
compactifying coordinate near infinity due to infinite blueshift in spatial frequency. 

Mapping infinity to a finite coordinate seems to require infinite resolution. However, a suit- 
able time transformation discussed in the next section provides a clean solution to this problem. 
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2.2. Hyperboloidal compactification 

The idea is to transform the time coordinate as in We introduce 

T=t -( x+ ih)> (9) 

With the compactification |5]l we get the Jacobian 

d T =d t , d x = (- 1 + C Q 2 )<9 T + n 2 d p , 
where we define Q := 1 - p. The advection equation in the new coordinates (p, t) reads 

d T u H — d„u — . 
C p 

This equation has the same form, up to an additional free parameter, as the advection equation 
in the original coordinates but the meaning of the coordinates is different. Solutions to the 
above equation in the bounded domain p e [0, 1] correspond to solutions to the original advection 
equation in the unbounded domain x e [0, oo). The free parameter C expresses the freedom to 
prescribe the characteristic speed in the compactifying coordinate and the number of cycles in an 
infinite domain. To see this, we write the solution (|7]i in the new coordinates 

u(p, T ) = - sin (2n(C Q + r)) . (10) 

The solution is depicted in figure [3] at t = for two values of C. The number of cycles on 
the domain is finite and depends on C. The wave is resolved evenly through the compactified 
domain. In such representations of the solution it should be kept in mind that lines of constant r 
do not correspond to lines of constant t. 




Figure 3: The solution at time r = 1 is plotted for two values of C. On the left panel we have C = 1 and on the 
right panel C = 5. The number of oscillations, and therefore the wavelength of the solution, can be influenced by the free 
parameter C. 



The idea to introduce a coordinate transformation of time in combination with compactifi- 
cation comes from general relativity [5]. The time function Q approaches characteristics of 
the advection equation asymptotically. In general relativity, infinity along characteristic direc- 
tions is called null infinity. Time functions whose level sets approach null infinity are called 
hyperboloidal because their asymptotic behavior is similar to the asymptotic behavior of stan- 
dard hyperboloids 013. To see this, consider the rectangular hyperbola on the (x,f) plane, 
t 2 - x 2 = C 2 , with a free parameter C. Shifting the hyperbola along the t direction by r gives 
(f - t) 2 — x 2 — C 2 . Introducing t as the new time coordinate we write 

T = t- VC 2 + x 2 . (11) 
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We plot in figure [4] two families of hyperbolae with C — 1 and C — 5. The asymptotes of these 
hyperbolae on x > are the characteristics of the advection equation, t = t—x, for any value of C. 
Time surfaces Q and (|9} share the same asymptotic behavior, hence the name hyperboloidal for 
such surfaces. A suitable compactification along hyperboloidal surfaces sets the coordinate loca- 
tion of null infinity to a time independent value (see [8 1 for a discussion of conformal and causal 
properties of hyperboloidal surfaces and compactification on asymptotically flat spacetimes). 




Figure 4: Rectangular hyperbola jl 1 1 for C = 1 and C = 5 depicted on the (x, f) plane. The hyperbolae are flatter near 
the origin for larger C but their asymptotic behavior in leading order is independent of C. 

It is useful to summarize the hyperboloidal compactification technique with general expres- 
sions. We introduce coordinates t and p via 

T = t-h(x), ( 12 ) 

The height function, h, must satisfy \dh/dx\ < 1 in units in which the speed of light (advection 
speed in the above example) is unity, so that r is a time function. We also require that the gradient 
of the function D. = Q(p) is nonvanishing at its zero set. The zero set of Q. corresponds to infinity 
in x. The coordinate transformations have the Jacobian 

Q? dh dQ. 

d, = d T , d x = -Hd T +—d„, where H := —(p), L:=Q-p— . (13) 

L ax dp 

The advection equation becomes in this notation 

d - u + n TK7 d P u = °' < 14 > 

(1 - H)L 



The specific choices |5]l and |9| satisfy 

Q 2 _ 1 
(l-H)L ~ C 



(15) 



In general, the time transformation must be chosen such that we have asymptotically in x, or 
equivalently as Q approaches zero 

1 - H ~ 0{Q}). (16) 

This relation formulates that the asymptotic boost of time surfaces needs to approach the speed 
of light as fast as the infinite compression of space to have a uniform outgoing characteristic 
speed in a compactified domain. This intuitive explanation suggests the names boost function for 
H, and compress function for £1 
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2.3. Wave equation 

The advection equation discussed in the previous section is a special example because its 
characteristics propagate in only one direction. A more representative example with incoming 
and outgoing characteristics is the wave equation, 

dju - d\u = 0. (17) 

The characteristics of fP7"| ) on a bounded domain are plotted on the left panel of figure [5] We are 
interested in the problem on the unbounded domain x e (— oo, oo). With ( p3] > we get 



Q2 I q2 /q2\ 

d * + [ 2HdTdp " T dp + dp{H)dr " dp \T) dp 



u = 0. (18) 



The transformation |9) is not the right choice for the solution domain. Instead, we choose the 
height function of standard hyperboloids, h{x) = VS 2 + x 2 . We prescribe an arbitrary coordinate 
location for null infinity by modifying the zero set of the compress function. We set 



II p 2 \ T II p 2 \ , TT 2Sp 



Q=- => L=- 1 + ^ and H= „ „ . (19) 

2\ S 2 J 2\ S 2 ) S 2 +p 2 V ; 

The parameter S determines the coordinate location of null infinity on the numerical grid. 
The unbounded domain x e (-oo, oo) corresponds in the compactifying coordinate top e (-5, S). 
The wave equation becomes 



2p 2 2 2S Q (3S 2 +p 2 )pD. 

3 T U + —d T d p U -a d p U+ s2+p2 d T U + S2 ( s2 +p 2) P " = ' (20) 



The equation evaluated at infinity, {p = +S }, takes the form 

d T [d T ± 2 d p ) u = 0, 

reflecting that both boundaries are outflow boundaries and do not require boundary conditions. 

The characteristics of the wave equation |20} are plotted on the right panel of figure [5] No 
boundary conditions are needed because no characteristics enter the simulation domain. Further- 
more, the outgoing characteristics leave the domain smoothly through the outflow boundaries. 

The initial value problem for the wave equation (jT7J» is related to the initial value problem 
for ( |T~8] > by time evolution as indicated in figure |4] This may be an undesired complication in 
practical applications. If we are interested in the evolution of compactly supported data, we 
can keep the t surfaces in an interior domain which includes the data and apply hyperboloidal 



compactification only in the exterior domain as discussed in section 2.5 

The above calculation is made for the one dimensional, source free, constant coefficient 
wave equation. The method applies when variable coefficients or suitable lower order terms 
are present. The requirement on the variable coefficients in the principal part is that they are 
asymptotically constant. For the lower order terms we require a fall off behavior of at least Q 2 
so that the division by 1 - H 2 in ( fT8| leads to a regular equation. 
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Figure 5: On the left panel we plot the characteristics for the standard wave equation ]17) on the bounded domain 
x £ [-10, 10]. On the right panel hyperboloidal compactification has been applied with infinity located atp = S = ±10. 
The standard wave equation on a bounded domain requires boundary conditions for the incoming characteristics from 
both boundaries. There are no incoming characteristics with hyperboloidal compactification. 



2.4. Hyperbolic systems 

Consider the linear, homogeneous system of partial differential equations with variable coef- 
ficients 

d t u = Ad x u, (21) 

where u = (u\(x, t), U2(x, f), . . . , u n (x, t)) T , and A is an n x n matrix that may depend on x. The 
transformation ( fl2| > with Jacobian ( [T3] l leads to 

Q 2 

(1 + //A)d T u = — A<9 p u. (22) 

Assuming that the time transformation has been chosen to satisfy ( fT6) >, we require that the poly- 
nomial remainder of det(l + HA) by 1 - H vanishes asymptotically. This is a condition on the 
asymptotic form of the elements of A. For example, taking n — 2 we write 

a 2\ a 22 

The asymptotic condition for the applicability of hyperboloidal compactification reads 

1 + an + #22 - 0i2«2i + 011022 = 0. (23) 

A typical example is the wave equatio n §Tl\ written as a first order symmetric hyperbolic 
system. The wave equation takes the form (|21|i in the auxiliary variables v = d,u and w - d x ii, 
with 

v\ . / 1 , 



w r 10 



The condition (23 i is satisfied. The transformed system reads 

^ u= (i^W("f -h) 8 ^ (25) 

Note that this equation is not the first order symmetric hyperbolic form of the transformed wave 
equation ([T8|. The particular choice (fl9]> leads to the regular system 



. 1 / -2Sp S 2 +p 2 \ 



The outer boundaries at p — ±S are pure outflow boundaries. 

As a further example consider the one dimensional Maxwell equations for the electric and 
magnetic fields (E, H) 

d,E = —d x H, d,H = —d x E, 

The electric permittivity e and the magnetic permeability fi may be point-dependent. The equa- 
tions have the form ( f2"Tj ) with 

l\ and A= ~^(° 

We get after hyperboloidal compactification 

= " ~ Wu- (26) 



{ep.-H 2 )L\ e H 

In vacuum outside a compact domain we have e = eo and // = //o, where eo and //o are the 
electric and the magnetic constants. We need to choose the asymptotic behavior of H such that 



yeo/A) — H ~ 0(Q. 2 ). Then the Maxwell equations behave similarly to the wave equation (125 i 
near null infinity. 

The example of Maxwell equations suggests that including lower order terms or variable 
characteristic speeds in a compact domain are straightforward in the hyperboloidal method as 
long as the asymptotic form of the equations are suitable. The asymptotic characteristic speeds 
need to be constant and lower order terms need to have compact support or fall off sufficiently 
fast. Hyperboloidal compactification can then be applied outside a compact domain as discussed 
in the next section. 



2.5. Hyperboloidal layers 

It may be desirable to employ specific coordinates in a compact domain without the time 
transformation or the compactification required by the hyperboloidal method. One reason is 
technical. Elaborate numerical techniques to deal with shocks, scatterers, and media assume 
predominantly specific coordinates. It may be impractical to modify these methods to work with 
hyperboloidal compactification throughout the simulation domain. Another reason is initial data. 
One may be interested in the evolution of certain (compactly supported) initial data prescribed 
on a level set of t. Therefore, it may be favorable to restrict the hyperboloidal compactification 
to a layer. 

We discuss briefly the perfectly matched layer (PML) by Berenger (9) to set the stage for 
hyperboloidal layers. In the PML method one attaches an absorbing medium — a layer — to the 
domain of interest such that the interface between the interior domain and the exterior medium 
is transparent independent from the frequency and the angle of incidence of the outgoing wave. 
Inside the layer the solution decays exponentially in the direction perpendicular to the interface. 
As a consequence, the solution is close to zero at the outer boundary of the layer where any stable 
boundary condition may be applied. The reflections from the outer boundary may be ignored if 
the layer is sufficiently wide. 

The success of the PML method lies in the transparency of the interface between the interior 
domain and the layer. This property finds explanation in the interpretation of Chew and Weedon 
of the PML as the analytic continuation of the equations into complex coordinates ifTUl . The 
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challenge is then to find suitable choices of the equations and the free parameters that lead to 
exponential damping of the solution in a stable way, which may be difficult depending on the 
problem Ol OH ED- 

A simple example omitting details of the method beyond our needs demonstrates the basic 
idea. We perform an analytic continuation of x into the complex plane beyond a certain interface 
R. The coordinate x is written in terms of its real and imaginary parts as Re(x) + icrlm(x)&(x-R), 
where <x is a positive parameter, and denotes the Heaviside step function. Setting k = 1, the 
plane wave ([!]) at time t — becomes 



u(x, 0) = e 



2mRe(x) -2wa-lm(x)@(x-R) 



(27) 



The strength of the exponential decay is controlled by a free parameter <x. The solution is plotted 
in figure|6]for <x = 0.1. 




Figure 6: The plane wave solution with the spatial coordinate analytically extended into the complex plane \27\ demon- 
strates how outgoing waves are damped exponentially in the absorbing layer of the PML method. We set cr = 0. 1 and 
R = 5. The dashed line indicates the location of the interface. 



For a hyperboloidal layer we perform a real coordinate transformation, both of space and 
time, beyond a certain timelike surface x — R called the interface. We set 



x-R- , £2=1 t®(p - R) . 



(28) 



where the coordinate location of infinity satisfies S > R. The width of the layer is S - R. This 
transformation is partly motivated by the agreement between the coordinates x and p along the 
interface in the sense that 



dx d x 

x(R)=R, ^-(R) = h -r^ 
dp dp L 



(R) = 0. 



(29) 



A simple prescription for the boost function is to require unit outgoing characteristic speed 
across the layer. For example, the outgoing and incoming characteristic speeds c ± for the wave 
equation on p > are 



Q 2 



with 



L = n-(p-R)- 



dQ, 



L{±\-H) r dp 

The requirement of unit outgoing characteristic speed reads c+ = 1 , implying 

O 2 

1 -H = — . 



(30) 



(31) 
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For p < -R we require c_ = -1. The resulting characteristics are depicted in figure [7] For 
|p| < R we obtain the standard characteristics in (x, t) coordinates. For \p\ > R we obtain the 
hyperboloidal characteristics in (p, t) coordinates (compare figure|5]l. 



T 




Figure 7: Characteristic structure for hyperboloidal layers with boundaries located at p = ±5 = ±10 and interfaces at 
p = ±R = ±5. Compare the inner domain [-5, 5] to the left panel of figure[5] and the layers [-10, -5] and [5, 10] to the 
right panel of figure]?] No characteristics enter the computational domain. 

The plane wave in the layer has the same form as in the interior by ( |28] l and ( [3Tj ) (figure [8] 
the dashed line indicates the interface). The important difference to the standard method is that 
the incoming characteristic speed vanishes at the outer grid boundary. 
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Figure 8: On the left panel we plot the plane wave solution with the compress function given in J28J and the boost 
function given in (31) . The characteristic speeds are depicted on the right panel in which the top and the bottom curves 
correspond to outgoing and incoming characteristic speeds. The domain is given by p e [0, 10] with a compactifying 
layer starting at p = R = 5 as indicated by the dashed lines. The incoming characteristic speed at the outer boundary 
vanishes, therefore no outer boundary conditions are needed in the layer. 



There is a large freedom in the choice of compress and boost functions, which may be ex- 
ploited for specific purposes. As an example, take a linear compress function 

0= l®(p-R), (32) 

j — K 

and the boost function of standard hyperboloids translated by R 

H - ——====®(x - R) . (33) 
V(* - R) 2 + C 2 

The resulting representation of the plane wave is plotted on the left panel of figure|9]for C — 3. 
Now we have fewer oscillations in the layer than in the interior because of the spatial redshift 
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controlled by C. A strong redshift, and consequently few spatial oscillations, may be preferable 
in a high order space discretization scheme in which a few grid points per wavelength are suf- 
ficient for good accuracy. On the right panel of figure |9| however, we see that strong spatial 
redshift comes at a price: the outgoing coordinate speed increases strongly in the layer. Evalua- 
tion of c+ from ([30} gives at infinity c+ = 2(S - R) 2 /C 2 . A small value for C with a wide layer 
leads to a high outgoing characteristic speed, which requires small time steps in an explicit time 
integration algorithm. 



ii Speeds 
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Figure 9: On the left we plot the plane wave solution with (32) and (33) for C = 3 in the otherwise same setting as in 
figure [8] There are fewer oscillations in the layer due to the stronger redshift causing a higher outgoing characteristic 
speed as depicted on the right panel. 

The balance between accuracy in time and in space is influenced by the compress and boost 
functions, and the free parameters included in them. The specific choices will depend on the 
requirements of the problem. It can be expected that in most applications the outgoing charac- 
teristic speed in the layer will be chosen close to the speed in the interior domain as in figure 

E 

The hyperboloidal layer is similar to the PML in that both methods allow us to solve the 
equations of interest in standard coordinates in an inner domain. In both methods the interface 
between the inner domain and the layer is transparent on the analytical level, independent of 
the frequency and angle of incidence of the outgoing wave. The key difference is that the PML 
absorbs the outgoing wave so that it is damped exponentially, whereas the compactifying layer 
transports it to infinity. The solution in the hyperboloidal layer is of interest, as opposed to 
the solution in the absorbing layer. In fact, the solution at the outer boundary of the layer is 
of special interest in radiative systems because it gives the radiation signal as measured by an 
idealized observer at infinity. 

3. Multiple dimensions 

Hyperboloidal compactification is applicable in multiple dimensions when compactification 
is performed in the outgoing direction. The main difference in multiple dimensions is that a 
rescaling of the unknown fields needs to be performed such that the fields attain a non-vanishing 
finite limit at infinity. The rescaling depends on the fall off behavior of these fields, and therefore 
on the equation and the spatial dimension. 
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3.1. Reseating 

We begin our discussion with the three dimensional wave equation on flat spacetime 



(-5? + A R3 )w = 0. 



Here, Ajy is the Laplace operator on three dimensional Euclidean space. We write the wave 
equation in spherical coordinates to single out the outgoing direction 

, , 2 1 \ 
-d, + d; + -d r + -rA S 2 m = 0, 
r r l I 

with A§2, the Laplace-Beltrami operator on the two sphere, and r, the radial coordinate. In sec- 
tion 2.3 we showed that hyperboloidal compactification introduces a divisor that is proportional 
to the square of the compress function. With the compactification r — p/Q we see that any term 
beyond the flat wave operator on the (r, t) plane that falls off as r~ 2 or faster is multiplied with at 
least Q 2 , and is therefore regular under hyperboloidal compactification. 

The angular part of the wave equation in spherical coordinates admits a regular compacti- 
fication. The first radial derivative term, however, results in a singular operator at infinity. We 
resolve this problem by rescaling. The three dimensional wave equation on the (r, t) plane takes 
the form of the one dimensional wave equation in the rescaled variable v := ru. We get 

-3 2 +3 2 + ^A s2 Jy = 0, 

which admits a regular hyperboloidal compactification. 

This procedure generalizes to dimensions other than three. The essential feature of the rescal- 
ing is that it takes care of the fall off behavior of the scalar field. In three dimensions ru attains 
a regular limit at infinity; in n dimensions the rescaling depends on n. Solutions to the wave 
equation decay asymptotically as r - ^ -1 ^ 2 by energy conservation. We expect therefore that the 
n dimensional wave equation written for the rescaled variable v := r" _1 ^ 2 M admits a regular 
hyperboloidal compactification. 

The n dimensional wave equation in spherical coordinates reads 



(-d 2 +d 2 r + "-^-dr + ^A s „-, J u = 0, 



where Ag»-i is the Laplace-Beltrami operator on the n — 1 sphere. Transforming to the rescaled 
variable v := r ( " _1 ^ 2 M, we get 

-d 2 +d 2 r - i(n - m - 3) + ^A s „-,Jv = 0. 

All terms beyond the one dimensional wave operator fall off as r~ 2 and are therefore amenable 
to a regular hyperboloidal compactification. 

The inclusion of sources or of suitable nonlinearities is straightforward under the condition 
that the corresponding terms in the equation fall off sufficiently fast. For example, a power 
nonlinearity on the right hand side of the wave equation of the type u p leads to the source term 
v P r -(«-i)(p-i)/2 -pjjjg term j s regular under hyperboloidal compactification if the power of r is -2 
or lower, which implies p > 1 + 4/(n - 1). The critical power for which equality is satisfied is 
also the critical conformal power for semilinear wave equations with a power nonlinearity. This 
apparent coincidence is explained in the next section. 
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3.2. Conformal method 

Compactification of spacetimes with a suitable time transformation as proposed by Penrose 
in [ 5 1 as well the hyperboloidal initial value problem as proposed by Friedrich [ 7 1 employ confor- 
mal methods. The conformal language is prevalent in studies of spacetimes in general relativity 
Ifl4l . In this section we discuss the hyperboloidal compactification using conformal techniques. 
This viewpoint is of theoretical and practical interest because it reveals the interplay of confor- 
mal geometry, partial differential equations, and numerical methods within the hyperboloidal 
approach, and also simplifies the implementation of the method for certain problems. 

In general, a hyperboloidal time transformation with a spatial compactification leads to a 
singular metric. Consider the Minkowski metric in spherical coordinates 

r, = -dt 2 + r 2 dr 2 + r 2 d(r 2 , (34) 

where da 2 is the standard metric on the unit sphere. Introducing new coordinates t and p as in 
(L2) gives 

o 2.HL 1 — H 2 2 2 2 

T] - -dr — -^-drdp + — — ^ — L dp + ^dcr . 

This representation of the Minkowski metric is singular at infinity. The singularity is removed 
by conformally rescaling the metric, 

1 — H 2 

g = Q 2 r] = -Q 2 dr 2 - 2HL drdp + —. L 2 dp 2 + p 2 do 2 . (35) 

Ll z 

The conformal metric g is regular at null infinity by ( fT6| ), and can be extended beyond null infinity 
in a process referred to as conformal extension ll5l fT4l[T5ll . In this context, the function Q is called 
the conformal factor. The zero set of the conformal factor corresponds to null infinity where it 
has a non-vanishing gradient. These properties of the conformal factor underlie our choices for 
the compress function in the previous sections. 

Partial differential equations within the conformal framework have first been studied for 
fields with vanishing rest mass, such as scalar, electromagnetic, and gravitational fields fl5l . 
The key observation is that the complicated asymptotic analysis of solutions to these partial dif- 
ferential equations is replaced with local differential geometry by considering the conformally 
transformed equations in a conformally extended, regular spacetime ifTBirPTl . 

We discuss the wave equation on Minkowski spacetime as an example. Under a conformal 
rescaling of the Minkowski metric, g — Q 2 77, the wave equation transforms as lfT4l[T5ll 

□„- — -R[g])v = Q.- ( " +3) ' 2 n n u, with v Q (1 -" )/2 M . (36) 
An j 

Here, O g := g' JV V A( V v is the d'Alembert operator with respect to g, R[g] is the Ricci scalar of 
g, and n is the spatial dimension of the spacetime. The rescaling with Q in the definition of the 
variable v is asymptotically equivalent to the rescaling in Section |3.1| where we factor out the 
fall off behavior of u such that the rescaled variable v has a non-vanishing limit at null infinity. 
To see this, consider a specific choice for O from previous sections, say Q(p) = 1 - p as in (|5). 
In terms of the coordinate r - p/O, the conformal factor reads Q(r) = (1 + r)~', which behaves 
asymptotically as r~'. Therefore the definition of v in < [36] > corresponds asymptotically to the 
definition of v in section [37X1 
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We can also explain the observation made at the end of section 3. 1 concerning the agreement 
between the critical conformal power and the critical power for which hyperboloidal compact- 
ification leads to a regular equation. Using D n u = u p and the definition of v in ([36) we get at 
the right hand side of the conformally invariant wave equation Q(("~ 1 )' : '~(" +3 W/ 2 v/> . Regularity of 
this term at infinity, where the conformal factor vanishes, requires (n - l)p > n + 3, which is the 
condition of Section [3~T] Equality is obtained for p c = 1 + 4/(n — 1) for which the semilinear 
wave equation is conformally invariant, hence p c is called the critical conformal power. 

The conformal approach may be useful for various reasons. It extends directly to asymp- 
totically flat backgrounds with non- vanishing curvature EJUS). It also helps identifying the 
transformation behavior of the equations independent of coordinates. For example, Yang-Mills 
and Maxwell equations are conformally invariant, and therefore do not require a rescaling of the 
variables (the specific variables in which the covariant equations are written may not be confor- 
mally invariant and may require a rescaling). 

Hyperboloidal compactification with rescaling, as presented in section 3.1 seems numer- 
ically feasible in any spatial dimension. In even spatial dimensions, however, the notion of 
conformal infinity may not be feasible fT9l . It is an open question whether the hyperboloidal 
technique applies to this case. If the method fails in even dimensions, it should be interesting to 
understand whether this failure is related to the violation of Huygens' principle. 

3.3. Hyperboloidal layers in multiple dimensions 

In multiple dimensional problems, the layer technique is employed along the outgoing direc- 
tion r. We employ a compactifying coordinate p defined via r — p/O. The compress function 
needs to be unity in a compact domain bounded by radius R, be sufficiently smooth across the 
interface at r = p — R, and vanish at a coordinate location S > R with non-vanishing gradient. A 
suitable choice is 

\S-RJ F ' (S-Kf F 

The coordinates r and p coincide up to second order along the interface (compare (|29|). The 
height function in the new time coordinate depends only on the radial coordinate: r = t - h(r). 
The boost function is given by H = dh/dr. It may be set such that the outgoing characteristic 
speed through the layer is unity. For the wave equation this corresponds to 

Q? 

1 -H = — , 
L 

as in pi) . Finally, the rescaling of the variable is performed with the c ompr ess function noting 



that asymptotically Q. ~ r . Further details are presented in section 4.2 in which the layer 
technique is applied to solve the wave equation numerically in three spatial dimensions. 

The extension of the hyperboloidal method from one dimension to multiple dimensions is 
fairly straightforward in spherical coordinates, which may be, however, too restrictive for the 
grid geometry in applications. The discussion of the conformal regularity of the Minkowski 
metric suggests that the method can also be applied in nonspherical coordinates as discussed in 
the next section. 
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3.4. Nonspherical coordinate systems 

Spherical coordinates are suitable for hyperboloidal compactification because cuts of null 
infinity have spherical topology fl5l . The coordinate shape of null infinity, however, does not 
have to be a sphere. We can employ any coordinate system with a unique outgoing direction for 
the compactification. 

An example for a nonspherical coordinate system, useful especially in electromagnetism, is 
the prolate spheroidal coordinate system. The relation between Cartesian coordinates {x, y, z] and 
prolate spheroidal coordinates {p, v, <f>} reads 

x = sinh// sin v cos^p, y — sinh// sin v sin <p, z- cosh //cos v. 

We have r 2 = sinh 2 p + cos 2 v. Here, p is the outgoing direction that has closed coordinate 
surfaces. Compactification is performed along p. Using the conformal method, we argue that 
if conformal compactification of Minkowski spacetime in these coordinates leads to a regular 
metric, the equations we solve on that background will be regular. The Minkowski metric reads 

77 = -dt 2 + dx 2 + dy 2 + dz 2 = -dt 2 + (sinh 2 /< + sin 2 v) (dp 2 + dv 2 ) + sinh 2 // sin 2 vdip 2 . 

We introduce a new time coordinate by setting 



t — t — Jl + sinh 2 //. 

The metric becomes 

rj - -dr 2 - 2 sinh p dpdr + sin 2 v dp 2 + (sinh 2 p + sin 2 v)dv 2 + sinh 2 p sin 2 vdip 2 . 
Compactification along p is performed via 

The conformal metric g - Q. 2 j] becomes 

g = -£l 2 dT 2 - Ipdpdr + sin 2 v dp 2 + (p 2 + Q 2 sin 2 v)dv 2 + p 2 sin 2 vdip 2 . 

The qualitative behavior of this metric near infinity is similar to the conformal Minkowski metric 
given in ( |35j ); the only difference is the coordinate representation. Therefore we conclude that 
hyperboloidal compactification of suitable hyperbolic equations in prolate spheroidal coordinates 
leads to regular equations as for spherical or oblate spheroidal coordinates. 

It is an open question whether the hyperboloidal method applies to Cartesian or cylindrical 
coordinates. The difficulty is to ensure regularity of the equations at corners and along edges 
with respect to limits to infinity. This requirement implies a restriction on the geometry of both 
the interface boundary and the numerical outer boundary. Even with the requirement of smooth 
coordinate surfaces, it would be useful to extend the method such that the layer has an arbitrary 
shape in coordinate space. This generalization would increase the efficiency of the method when 
dealing with scatterers with irregular shape. 
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4. Numerical experiments 



The aim of this chapter is to demonstrate numerical implementations of hyperboloidal com- 
pactification on simple examples in one and three spatial dimensions. The application of the 



method in two dimensions is an open question as discussed in 3.2 and is left for future work. 
Also left for future work is an in-depth numerical analysis of the method in comparison with 
other boundary treatments, such as absorbing boundary conditions or perfectly matched layers. 

In previous sections, we have seen that we may employ an arbitrary coordinate system in an 
interior domain restricting compactification to a layer outside that domain. This idea is based 
on the matching method with a transition zone presented in |S). Numerical applications of the 
matching, however, are inefficient due to a blueshift in frequency in the matching region and 
many arbitrary parameters for the transition function [20, 2T|. A hyperboloidal layer gives a 
sufficiently smooth interface without a transition function. 

In this section, we compare the numerical accuracy of solutions with and without the layer 
in one dimension. In three dimensions we focus on the application of hyperboloidal layers for 
the wave equation. Calculations using hyperboloid foliations, that is, constant mean curvature 
surfaces without the layer, have been presented in ll2"Tll . 



4.1. One spatial dimension 
4.1.1. Analytical setup 

Consider the Maxwell equations ( p6*| i. Assume that the electric permittivity and the magnetic 
permeability are constant and have unit value. The transformed Maxwell system for the unknown 
vector u = (E, ft) reads 

/ H 1 

(l-H 2 )L\ 1 H 



drn = ~ „ rrtw | , „ \d p u. (37) 



The characteristic speeds are c± = -Q 2 /((+ 1 + H)L). This system is similar to the wave equation 
written in first order symmetric hyperbolic form ( |25) , so our results apply both to Maxwell and 
wave equations in one dimension. 

We experiment with two sets of the compress and boost functions. First we employ the 
hyperboloid foliation everywhere in the simulation domain. We set as in ([19} 

a=Ul-£\ and H= 2SP 



2\ S 2 j S 2 +p 2 ' 

The characteristic speeds read c ± = +(1 + p/S) 2 /2. The minus sign corresponds to the incoming 
speed at the right boundary, which vanishes at p = S . The plus sign corresponds to the incoming 
speed at the left boundary, which vanishes at p = -S . The qualitative behavior of the charac- 
teristics is the same as depicted on the right panel of figure [5] The time step in an explicit time 
integration scheme is restricted by the maximum absolute value of the characteristic speed which 
reads |c| max = 2. 

The second set of compress and boost functions are for hyperboloidal layers. We set the 
compress function as in ( |2"8j ) 

(p - R) 2 
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We determine the boost function from the requirement of unit outgoing characteristic speeds 
through the layers. We set 

Q. 2 Q} 

H = 1 for p > R , and H = -1 for p<-R. 

L L 

The characteristic structure for the resulting equations are depicted in figure [7] 

4.1.2. Numerical setup 

The hyperboloidal method is essentially independent of numerical details due to its geometric 
origin. We discretize ([37) employing common methods: an explicit 4th order Runge-Kutta time 
integrator and finite differencing in space with 4th, 6th, and 8th order accurate centered operators. 
At the boundaries we apply one-sided stencils of the same order as the inner operator. We use 
Kreiss-Oliger type artificial dissipation to suppress numerical high-frequency waves [22]. For a 
2p — 2 accurate scheme we choose a dissipation operator of order 2p as 

h 2 P- 1 

Ai lss = ei-lY^D^Di, 

where h is the grid size, D ± are forward and backward finite differencing operators and e is the 
dissipation parameter. 

Both for the hyperboloid foliation and the layer we set S = 10. The simulation domain is 
then given by p e [-10, 10], which corresponds to the unbounded domain x e (-00,00). The 
interface for the layer is at R = ±5. Hyperboloidal coordinates (p, t) coincide with standard 
coordinates (x, t) within the domain |p| < 5. 

We solve the initial value problem for ( |3"7| i with a Gaussian wave packet centered at the origin 
for the electric field and vanishing data for the magnetic field. We set at the initial time surface 

E(p,0) = e- p2 , fl(p,0) = 0. 

The solutions with respect to the hyperboloid foliation and the hyperboloidal layer correspond 
to different initial value problems due to the different time surfaces. The solution constructed 
using the hyperboloidal layer corresponds, however, to the solution that we would obtain using 
the standard coordinates (x, t). 



4.1.3. Results 



Figure 10 shows convergence factors for the electric field from a three level convergence test 
with 100, 200, and 400 grid cells. The convergence factor Q is calculated for the electric field in 

\\ plow _ pined \\ 

the Li norm via Q := log 2 p^arpaj ■ The factors are in accordance with the implemented finite 
difference operators. 

The total energy is radiated to infinity leaving the zero solution behind in accordance with 
Huygens' principle for Maxwell equations in one dimensional fiat spacetime. Therefore, a good 
measure of the quality of the boundary treatment is the value of the unknown after the initial wave 
packet leaves the simulation domain. This value is related to the numerical reflection coefficient 
of the boundary. The analytical reflection coefficients at the interfaces and at the boundaries are 
zero in the hyperboloidal method. Numerically, however, the interfaces and the finite differencing 
at the outer boundaries cause reflections. Consider the L2 norm of the solution as a function of 
time plotted in figure [TT] The plots depict the field with respect to the hyperboloid foliation (left) 
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Figure 10: Convergence factors in time in the Lj norm for the Maxwell equations in one spatial dimension. Solid curves 
represent the hyperboloid foliation, dashed curves the layer method. The curves indicate from bottom to top 4th, 6th, and 
8th order convergence in accordance with the order of implemented finite differencing operators. 



and with respect to the layer (right). The errors are larger for the layer. Red, green, and blue 
curves represent solutions calculated with 4th, 6th, and 8th accurate finite difference operators, 
respectively. To each order the solution is calculated with and without dissipation. We observe 



that artificial dissipation strongly reduces numerical errors. Figure 11 shows that the error at late 



times with 4th order stencils with dissipation has the same order of magnitude as the one with 
8th order stencils without dissipation. 
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Figure 1 1 : The L2 norm of the electric field in time for 4th (red), 6th (green), and 8th (blue) order finite differencing with 
respect to the hyperboloid foliation (left) and the hyperboloidal layer (right). The curve indicating smaller errors to each 
order corresponds to the solution with artificial dissipation. The plot shows that the effect of dissipation to the quality of 
the solution is comparable to using a higher order numerical discretization. Errors with the layer method are larger than 
eiTors with the hyperboloid foliation. 



The main result from figures 10 and 11 is that the numerical error is reduced drastically by 
higher resolution, higher order discretization, and artificial dissipation — without changing the 
boundary treatment. The boundary treatment does not introduce errors into the solution that are 
put in by hand. To emphasize this point, we plot in figure 12 two solutions using the hyperboloid 
foliation and the hyperboloidal layer. For each method, one solution is obtained with 4th order fi- 
nite differencing, without dissipation, and 100 grid cells (red curve); the other solution is obtained 
with 8th order finite differencing, with dissipation, and 200 grid cells (blue curve). The combined 
effect of these improvements demonstrates that accuracy with hyperboloidal compactification is 
not restricted by the boundary treatment but by the numerical accuracy in the simulation domain. 
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Figure 12: Comparison of numerical errors as indicated by late time values of the Li norm of the electric field between 
a numerical solution with low (red) and high accuracy (blue) for the hyperboloid (left) and the layer (right) method. 
Accuracy with hyperboloidal compactification is not restricted by the boundary treatment. 



4.2. Three spatial dimensions 
4.2.1. Analytical setup 

Consider the semilinear wave equation with a focussing power nonlinearity, D n u = — u p . We 
solve this equation with and without the nonlinear term by using the conformal method. The 
semilinear conformal wave equation reads 

□ g v= ^R[g]v-Q?- 3 v p , (38) 

where g is the conformal metric, R[g] is the Ricci scalar to the metric g, and v = Q. l u. The 
conformal metric g is obtained from the Minkowski metric rj given in ((34]) by the time transfor- 
mation t — t — h(r), the spatial compactification r - p/O, and the conformal rescaling g = Q. 2 ?]. 
We set the conformal factor Q. as in section [331 



\S-R) ^ (S-R) 4 ^ 

This choice ensures that the coordinates r and p coincide and the Ricci scalar is continuous along 
the interface. The boost function is such that the outgoing characteristic speed through the layer 
is unity pi) . The conformal metric becomes 

g = -Cl 2 dT 2 - 2(L - Q 2 ) drdp + (2L - Q 2 ) dp 2 + p 2 dcr 2 . 

For p < R we recover the Minkowski metric ( |34) . The characteristic speeds of spherical wave 
fronts are 

L + L-D? 



2L-Q. 2 ' 

We observe that the outgoing characteristic speed c + is unity everywhere; the incoming charac- 
teristic speed c_ is unity inside the domain p < R, and zero at infinity. The characteristic structure 
is similar to the right half of figure [7] The Ricci scalar reads 

6Q(OL' - 21X1') 
R[8] p 2 l? ■ 

The apostrophe denotes derivative by the argument. 
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4.2.2. Numerical setup 

We apply similar numerical techniques as those that have been used to test constant mean 
curvature foliations Ell . The application of hyperboloidal compactification is complicated by 
the requirement of a spherical grid outer boundary because spherical coordinates are ill-defined 
at the origin. To deal with this problem, we use the Spectral Einstein Code (SpEC) l23ll . 




Figure 13: The numerical grid for calculations in Minkowski spacetime. We have a cube around the origin on the domain 
[-2, 2] in each Cartesian direction and four spherical shells extending from p = 2 to future null infinity at p = 20. The 
colors depict off-centered Gaussian initial data for the time derivative of the scalar field. 



We cover the origin with an inner cube using Cartesian coordinates, thereby avoiding the 
coordinate singularity of spherical coordinates at the origin. Spherical shells starting within the 



cube extend to a spherical outer boundary which corresponds to null infinity (figure 13 1. The 
domain in the radial direction is p e [0,20]. The cube around the origin has x, e [-2,2]. Four 
spherical shells extend from p = 2 to future null infinity. 

Spatial derivatives are discretized by a pseudospectral method. For the Cartesian grid we use 
Chebyshev polynomials in each direction. For the shells we use Chebyshev polynomials in the 
radial direction and spherical harmonics in the angular directions. The wave equation is written 
in first order symmetric hyperbolic form and time integration is performed with a Runge-Kutta 



algorithm (see ETll for details). The colors in figure 13 depict the off-centered data prescribed 
for the time derivative of the scalar field. The data for the scalar field vanishes. 

The interface to the hyperboloidal layer is at R - 10. Experiments confirm that the results 
reported below are not sensitive to the location of the interface if the layer is sufficiently thick to 
resolve an outgoing wavelength. The optimal thickness of the layer in applications will depend 
on certain factors such as the wavelength of the outgoing radiation or, probably, the shape of the 
scatterer. Future research should determine guidelines for the thickness of the layer. 

4.2.3. Results 

Spatial truncation errors converge exponentially in a pseudospectral code. We show such 



spectral convergence in figure 14 for a solution of the homogeneous linear wave equation with 



off-centered initial data. We calculate convergence in the auxiliary variable of the first order 



reduction v, := 3, v. Figure 14 shows that the numerical difference between v, and d, v converges in 
the L2 norm geometrically with N collocation points in each subdomain direction. The constraint 
error grows slowly in time but the evolution is stable. Convergence of the constraint error is not 
disturbed by the presence of the boundary located at p = 20. 
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Figure 14: Spectral convergence for the homogeneous wave equation with oft-centered initial data. We show the L2 norm 
of the constraint, Vj - d,v, using N = 5, 7, 9, 1 1 collocation points (from top to bottom) in each subdomain direction. 



A stringent test for boundary methods is the inclusion of nonlinear terms in the equations. 
The boundary treatment for semilinear wave equations is a difficult problem, especially in three 
spatial dimensions 0). Therefore, we present also results for the wave equation with a cubic 
source term (with p — 3 in (|3~8|l). The backscatter due to the self interaction of the field violates 
Huygens' principle which makes the solution more interesting. Transparency boundary condi- 
tions must not eliminate all reflections from the outer boundary but only spurious ones. This 
makes the treatment of the artificial boundary difficult. 

The backscatter by the nonlinear term leads to a late time polynomial decay of the solution: 
as r 2 at a finite distance, as r 1 at infinity. It is difficult to obtain the correct decay rates with an 
artificial outer boundary. The decay rate at infinity is not even accessible with standard methods. 



Figure 15 shows that the rates are calculated accurately with the hyperboloidal method. The 



local decay rate is defined as d\n\v(p, r)|/c/lnr. The invariance of the time direction under 
hyperboloidal transformations implies that the rates calculated with the hyperboloidal method 
are equivalent to the rates calculated with the untransformed equations. 
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Figure 15: Local decay rates for the cubic wave equation measured by observers located from top to bottom at p = 
(20, 19.97, 19.9, 19.8, 19.4, 17.861, or equivalently at r = {00, 1970, 500, 226, 87, 29}. The bottom curve that corresponds 
to the fastest decay bends up after t = 500 due to accumulation of numerical errors. 



Any numerical method fails after a certain time due to accumulation of numerical errors. 



In figure 15 the bottom curve that corresponds to the fastest decay bends up after about r = 
500. This behavior is expected, and is due to accumulated numerical errors. We can delay its 
appearance by increasing the numerical resolution of the simulation. The long time accuracy of 
the solution is not limited by the boundary treatment, but solely by the accuracy of the interior 
calculation. 
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5. Discussion 

Hyperboloidal compactification provides a clean solution to the artificial outer boundary 
problem. It allows us to solve suitable hyperbolic equations accurately on unbounded domains. 
The underlying idea originates in studies of asymptotic structure of spacetimes in general rela- 
tivity 13 . It is remarkable that such a practical method is available at the interface between 
differential geometry, conformal structure, partial differential equations, and numerical analysis. 

With the hyperboloidal layer method, in which compactification is restricted to a layer, we 
can apply arbitrary coordinates in an interior domain, and we gain direct quantitative access to the 
asymptotic solution. The numerical boundary of the simulation domain corresponds to infinity, 
therefore no outer boundary conditions are needed. The hyperboloidal time transformation leads 
to a non-vanishing coordinate speed for outgoing characteristics up to and including infinity in 
compactifying coordinates. 

Some further advantages of the hyperboloidal method are listed as follows. Efficiency of nu- 
merical calculations improves by suitable choices of free parameters that modify the wavelength 
of outgoing radiation. It is straightforward to apply the method to various covariant hyperbolic 
equations, such as wave, Maxwell, or Yang-Mills equations. It is not necessary to calculate 
boundary data that typically depend on the particular equations. Including source terms and 
nonlinearities does not modify the boundary treatment if certain asymptotic fall-off conditions 
are satisfied. No separate error controlling at the outer boundary is needed because there are no 
errors beyond the numerical discretization. There is no overhead in software implementation of 
boundary routines. The method is largely independent of numerical schemes due to its geometric 
origin. 

An important current limitation of the method is the requirement of a spherical or spheroidal 
grid near the outer boundary. It is not clear whether hyperboloidal compactification works with 
global Cartesian or cylindrical grid structures. Another limitation is source terms that lead to a 
singular compactification (Klein-Gordon equation, low power nonlinearities). The applicability 
of the method in even spatial dimensions, starting with two dimensional problems, should be 
studied. 

The decision to apply hyperboloidal compactification can be made only after detailed com- 
parative studies which go beyond the scope of this paper. A comparison of different boundary 
treatments (absorbing boundary conditions, perfectly matched layers, and hyperboloidal layers) 
within the same numerical setting on a simple example, such as the linear, homogeneous wave 
equation, would be useful. The method should be applied to Maxwell equations and its effi- 
ciency should be studied when scatterers with irregular shapes are present in the computational 
domain. Generalizations of the method, such as adapting the coordinate representation of the 
hyperboloidal layer to the shape of scatterers, might lead to very efficient computations. 

The relevance of hyperboloidal compactification depends on the problem at hand. The 
method in its current form does not apply to non-covariant problems, however, it may be ex- 
tended to deal with heterogeneous media, anisotropic elastic waves, optical waveguides, Euler 
equations, Navier-Stokes equations, among others. 

The idea of boosting the time direction may find applications beyond compactification. It 
seems that transformations of the characteristic cone for hyperbolic equations are largely unex- 
plored outside numerical relativity. The presented solution of the outer boundary problem may 
be just one application of such transformations. Others might be useful, for example in numerical 
studies of high frequency waves. 
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